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equation 


to  be  minimized  for  |u|  <_ m and  6=0.  By  considering  the  problem  as  € ->  0, 
singular  perturbation  concepts  can  be  used  to  compute  solutions  consisting  of 
bang-bang  controls  followed  by  singular  arcs.  The  procedure  further  develops  a 
numerical  technique  proposed  by  Jacobson,  Gershwin,  and  Lele,  as  well  as  ad- 
ditional analytic  methods  developed  by  other  authors. 


1.  INTRODUCTION 


A typical  singular  optimal  control  problem  consists  of  a state  equation 

j 


(1.1) 


x=Ax+bu,  0<t<T<” 


J 

i 


(subject  to  end  conditions  on  the  n - vector  x)  and  a scalar  cost  functional 


(1.2) 


uT 

= 1 VQx 
n 


which  is  to  be  minimized  for  a symmetric,  positive  semi-definite  matrix  Q 
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The  known  theory  of  singular  problems  has  recently  been  summarized  in  Bell  and 
Jacobson  (1975),  while  research  papers  on  the  general  problem  and  particular 
applications  abound  in  the  current  literature. 

Jacobson  and  coworkers  (cf.  Jacobson  and  Spe.yer  (1971)  and  Jacobson, 
Gershw'n,  and  Lele  (1970))  converted  such  problems  to  nonsingular  (though  nearly 
singular)  problems  by  considering  the  perturbed  cost 


0.4) 


J (£)  = \ (x'Qx  + £^u2)dt 
J0 


for  a sequence  of  positive  £ values  tending  toward  zero.  This  device  resulted 
in  considerable  progress  both  for  theoretical  purposes  and  for  numerical  calcula- 
tion. More  recently,  O'Malley  and  Jameson  (1975,  1976)  and  O’Malley  (1976)  have 
used  an  analogous  method  to  analytically  solve  such  singular  control  problems 
with  time-varying  coefficients  and  vector  controls,  but  without  bounds  on  the 
components  of  the  control.  Their  results  rely  heavily  on  the  asymptotic  theory 
for  singularly  perturbed  boundary  value  problems  for  ordinary  differential 
equations  (cf.  Wasow  (1965),  Lions  (1973),  and  O'Malley  (1974)).  Applications 
of  such  analyses  in  control  have  recently  been  surveyed  by  Kokotovic  et  al. 

(1976).  Bell  and  Jacobson  (1975)  stress  the  need  and  difficulty  of  developing 
computational  techniques  for  singular  control  problems  and  observe  that  the 
best  current  schemes  are  the  "£ -method"  of  Jacobson  et  al. (1970)  and  gradient 
techniques  (cf.,  e.g.,  Pagurek  and  Woodside  (1968),  Sirisena  (1974),  and  Edge 
and  Powers  (1976)).  We  note  that  these  £ -methods  are  analogous  to  the  arti- 
ficial viscosity  methods  popular  in  computational  fluid  dynamics  (cf.  Richtmyer 


and  Morton  (1967)).  In  addition,  we  note  that  Boggs  (1976)  has  overcome  some 
difficulties  with  gradient  techniques  by  using  asymptotic  methods. 
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Jacobson  et  al , (1970)  and  O'Malley  and  Jameson  (1975)  show  that  although 
the  optimal  cost  J(£)  given  in  (1.4)  converges  uniformly  as  £ -+■  0,  the 
limiting  control  generally  converges  nonuniformly.  Indeed,  the  control  must 
be  expected  to  consist  of  bang-bang  arcs  and  singular  arcs  (cf.  Johnson  and 
Gibson  (1963))  and  convergence  is  necessarily  nonuniform  as  £ -*■  0 at  switch 
points.  [In  another  context,  we  note  that  the  natural  sluggish  transfer  of 
mechanical  systems  at  switch  points  of  bang-bang  control  could  be  modelled  by 
another  singular  perturbation  device,  viz.  the  introduction  of  a cost  term 
'To*  2 

y (u(t))  dt  for  y small.]  We  also  recall  that  Powers  and  McDanell  (1971) 
and  Edgar  and  Lapidus  (1972)  report  practical  success  in  using  the  € -method  for 
Saturn  rocket  guidance  and  chemical  reactor  problems.  Practical  difficulties 
certainly  remain,  however.  As  € -»•  0,  for  example,  Jacobson  et  al  . (1970) 
note  that  numerical  instability  manifests  itself  and  they  suggest  a sufficiently 
good  approximation  might  result  by  reducing  £ to  a small,  but  "still  numerically 
stable"  value.  Analogous  difficulties  have,  of  course,  been  common  in  the 
numerical  solution  of  boundary  value  problems  for  stiff  ordinary  differential 
equations  (cf.  Willoughby  (1974)  or  Aziz  (1975)).  Substantial  progress  has 
been  made  on  these  problems  by  using  singular  perturbation  ideas  to  develop 
specially-tailored  numerical  methods  (cf.  Miranker  (1973),  Ferguson  (1975), 
Flaherty  and  O'Malley  (1977),  and  Kreiss  and  Nichols  (19  )).  We  propose  a 

similar  program  to  further  study  singular  control  problems.  By  combining 
asymptotic  and  numerical  ideas,  improved  methods  will  necessarily  follow. 

We  cannot,  of  course,  obliterate  the  unavoidable  complicated  behavior  inherent 
in  these  singular  problems  (cf.,e.g.,  Fuller  (1963)).  We  observe  that  the 
nearly  singular  problems  are  likely  to  be  of  independent  interest  (cf.  O'Malley 


| 
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and  Jameson  (1975)),  though  we  shall  not  pursue  that  question  here. 

In  their  recent  study  of  singular  arc  problems  with  unbounded  controls, 
Jameson  and  O'Malley  found  considerably  different  behavior  in  a sequence  of 
cases.  Case  k,  k = 1,  2,  ...»  corresponds  to  the  more  familiar  classification 
of  singular  arcs  of  order  k (cf.  Robbins  (1967)  or  Goh  (1966))  and  is  defined 
such  that  b ' ( A * )^QA^b  = 0 for  j = 0,  1,  ...,  k - 2 and  b'(A,)k~1QAk"1b  > 0. 

A given  problem  need  not  fit  into  any  such  case  (as  when  Q = 0)  and  for 
vector  controls  may  lie  between  cases  (when  the  last  matrix  is  singular,  but  of 
positive  rank  (cf.  Anderson  (1973)).  Like  much  of  the  literature  (cf.  Wonham 
and  Johnson  (1964),  Sirisena  (1968),  and  Anderson  (1972)),  our  discussion  will 
mostly  concern  Case  1.  For  Case  k problems,  the  singular  arc  solution  for 
an  n - vector  x and  scalar  u involves  a control  law  for  a dynamical  system 
of  order  n - k (cf.  O'Malley  and  Jameson  (1976)).  It  is  natural  then  to  seek 
a control  u which  switches  k - 1 times  between  control  bounds  +m  before 
reaching  the  singular  arc.  It  is  easy  to  see  that  such  a control  will  not 
generally  be  optimal  for  k >_  2,  though  we  show  that  it  is  quite  satisfactory 
for  certain  examples.  We  must  anticipate  such  difficulties,  however,  since 
experts  conjecture  that  infinite  switching  (as  in  the  classical  Fuller  problem) 
is  generic  for  k >_  2 (cf.  Krener  (1976)).  We  note  that  the  nonsingular 
problem  with  a small  fixed  value  of  £ could  have  an  optimal  solution  with  a 
finite  number  of  switchings,  while  the  limiting  singular  problem  involves  infinite 
switchings. 
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2.  TWO  FIRST  ORDER  SINGULAR  ARC  PROBLEMS 


a.  The  simplest  singular  arc  problem  may  be 


6 

'x  = u,  x(0)  = 1 
with 

1 (2  p 

0 = y x^(t)dt 
£>0 

^to  be  minimized  for  |u|  <_  m. 


For  m = ®,  the  optimal  control  features  an  initial  negative  delta  function 
impulse  which  drives  the  state  immediately  to  zero  with  the  optimal  cost  J*  = 0, 
To  solve  this  problem  by  the  6 -method  we  consider  the  nonsingular  problem 


' x = u,  x(0)  = 1 
wi  tti 


\ 


2 

J(£)  = y (X2(t)  + £^U^(t))dt 
C>  0 


to  be  minimized. 


Its  unique  solution  for  m =*  », 

jx(t.e)  - (i 

|u(t.e)  . 4(1  ♦ e-^rhe-^  - 
has  the  asymptotic  limit 

x(t,£)  % e and  u(t,t)  * - ^e 

for  € -*  0+.  We  note,  in  particular,  that  this  limiting  control  behaves  like 


It  f(t)e~t/edt  + f (0)  as  € -*  0. 
J0 


-6  at  t = 0 since 


For  bounded  controls,  |u | _<  m,  the  optimal  control  does  the  best  it  can  (cf. 
e.g.,  Sage  (1968)),  viz.  for  rn  > ~ 


■ 


The  integrated  effect  is  the  same,  i.e.,  the  singular  arc  solution  x = u = 0 
is  ultimately  reached  (see  Figure  1).  Indeed,  the  preceding  results  can  be 
recovered  by  letting  m ®.  For  bounded  m,  the  problem  can  be  explicitly 
solved  by  using  the  cost  functional  J(£)  and  seeking  the  asymptotic  solution 
as  £ -*■  0.  For  m = 1 , the  limiting  cost  is  exactly  ^ compared  to  Jacobson 
et  al  - (1 970) ' s calculated  values  0.1717  and  0.1617. 

Slight  modifications  of  this  example  can  be  easily  handled,  e.g.,  Rozonoer's 
example 

• 2 

x = u,  x(0)  = 1 , x(2)=p  | u | <_  1 , J = *2dt 


(cf.  Pagurek  and  Woodside  (1968))  and  the  vector  control  problem 


xi  = iu.,  |u.|  < 1,  i = 1,  2,  J = \ 


(x2  + x2)dt. 


b.  As  a second  example,  consider  the  harmonic  oscillator  problem 
f 

X-j  — X2  * Xr|(0)  — 0 

x2  = u,  x2(0)  = 1 
(2.1)  | with  cost 

J(€) 


^ (x2  + x2  + £2u2)dt 

0 1 *■ 


and  I u|  < m. 


For  m =»,  the  usual  state-costate  formulation  (cf.,e.g.,  Athans  and 
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Fa  1 b (1965))  implies  that  the  optimal  control  will  be  given  by 

(2.2)  u = -p2/e2 

where  the  state  and  costate  vectors  x = (x^,x2)'  and  p = (p^,p2)'  satisfy 
x i — x2,  x-|(0)  — 0,  x 2 — u j x2(0)  — 1 

(2.3) 

I - ~X-|  > P^(^)  ~ 0;  P2  - -x2  P-|  > P£(^)  ~ 

The  asymptotic  solution  of  the  singularly  perturbed  two  point  boundary  value 
problem  (2.2)-(2.3)  can  be  readily  obtained.  It  consists  of  an  initial 
boundary  layer  (endpoint  region  of  nonuniform  convergence)  while  the  limiting 
solution  within  (0,5)  satisfies  the  limiting  problem 

X10  = X20’  0 = p20’  p10  = ~X10’  and  p20  = ”X20  " p10 

obtained  when  we  set  £ = 0.  Since  p2q  = 0 and  X20  _p10’  we're  left  with  a 

linear  system  for  X-jg  and  p^g.  If  we  now  use  the  boundary  values  X-|g(0) 

= x-j(O)  = p ^ q ( 5 ) = P-|(5)  = 0,  we  get  the  trivial  singular  arc  solution 

(2.4)  X]0  = p10  = X2Q  = p20  = 0. 

Further  calculations  (cf.  O'Malley  and  Jameson  (1975))  also  show  that  the  limiting 
control  has  a negative  delta  function  impulse  at  t = 0.  We  note  that  the  limiting 
solution  for  t > 0 wouldn't  be  trivial  if  x-|(0)  1 0 and  that  introduction  of 
the  £ provides  a more  convenient  method  of  finding  the  limiting  singular  arc 
solution  than  the  more  familiar  technique  of  differentiating  the  optimality  con- 
dition Hu  = 0 twice  with  respect  to  t (cf.,e.g.,  Robbins  (1967)). 
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The  preceding  solution  wouldn't  be  appropriate  for  a finite  control  bound 
m,  because  the  impulse  in  the  initial  boundary  layer  would  exceed  this  bound 
Instead,  standard  maximum  principle  arguments  (cf.  Boltyanskii  (1971))  show 
that  the  continuous  optimal  control  is  determined  by 


'u  = - P2^^  if  [ P2 1 < m€^ 
and 

\ 

| u | = m otherwise. 

\ 

Moreover,  the  state  and  costate  vectors  satisfy  the  canonical  equations  (2.3) 
as  before.  Anticipating  that  the  control  initially  saturates  at  its  negative 
bound,  we  might  seek  a solution 


j 

I 

I 

f 

I 


* 


(2.5) 


u = 


-m. 


-p2/6‘ 


0 < t < ti 

i-}  £ t < 5 


and  determine  t-j , if  possible,  so  that  saturation  does  not  occur  on  (t-j  ,5) 
2 

while  | p2 1 > on  (0,t-|).  Clearly, 

(2.6)  x-|  ( t ) = - jt2  + t,  x2(t)  = -mt  + 1 

for  0 £ t < t-j  while  we  must  satisfy  the  singularly  perturbed  problem 


i it 

(2.7) 


VI 


/ • 

X-j  = X2,  X-j  ( t-j  ) = x i ( 1 1 ) 

2 • 

^ x2  - ~P2>  x2(t^ ) - X2(t-| ) 

< • 

Pi  x-| , Pi  ( 5 ) = 0 

’ 2 
P2  - _x2  ~ P] 5 P2  ( ) ~ ^ » 

\ 


p2(5)  = 0 


I 


'U.  S". 
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The  system  (2.7)  has  eigenvalues  +l/k(£),  +k(6)/c,  where  k2(6) 

= (1+/T  4£2)/2.  Thus,  the  solution  of  the  two-point  problem  is  of  the 
form 


(2.8) 


'x11 

-1/k  ' 

'1/k  ' 

x2 

-1/k2 

(t-t-,  )/k 

-1/k2 

= 

Ae  1 + 

pi 

1 

1 

% 

e2/k3 

-e2/k3 

Be 


)/k 


' -€/K2  ' 

-k(t-t,)/6 
Ce  1 + € 

r-  ,,  2\ 

-<z/k 

1/k 

-1/k 

2 3 

€?/k3 

e 

6 

j 

De-k(5-t)/6 


The  five  boundary  conditions  in  (2.7)  provide  four  linear  equations  for  A,  B, 

C,  and  D and  a nonlinear  equation  for  t-| . For  5 - t^  » 6,  the  exponential 
-k(5-t1 )/6 

e is  exponentially  negligible,  so  we  successively  find 

-2(5-t7)  -(5-t,) 

A -Be  1 , D « 2Be  1 , 

-2(5-t,)  , 

C « m + B(1  + e ) , and  B ^(x^t^)  - x2(t-|)) 


while  t-j  must  satisfy 

-2(5-t] ) mt2+2(m-l)t-,-2 

mt^-2(m+l  )t-|+2 


It  follows  that  there  is  a unique  root  t-j  such  that 


0 < t-l  < — 

i m 


- 1 + /7 


2‘ 


5EZ 


m 
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Indeed  for  m = 1,  is  very  near  the  upper  bound  /2”,  so  the  optimal 

control  is  approximately  given  by 


(2.9) 


0 < t < /2 


(/2  - 1 )e 


-(t-/2) 


t > /2. 


lie  solved  problem  (2.1)  with  6=0  by  our  asymptotic  technique,  which  is 
discussed  in  Section  3,  and  found  that  our  results  (see  Figure  2)  qualitatively 
agreed  with  the  computed  solution  pictured  in  Jacobson  et  al.  (1970).  We 
found  a minimum  cost  of  0.379  compared  to  their  value  of  0.414.  Anderson 
(1972)  solved  this  problem  by  a search  technique,  obtaining  a switching  time 
of  1.414  with  corresponding  trajectory  values  of  0.4144  and  -0.4136,  com- 
pared to  our  values  of  1.414,  0.4144,  and  -0.4138,  respectively. 

3.  A SYSTEMATIC  APPROACH  TO  CASE  ONE  PROBLEMS 
Let  us  now  consider  free  endpoint  problems  where 


(3.1)  b'Qb  > 0 

for  a prescribed  initial  state  vector  x(0).  Since  the  vector  b has  rank  one, 
we  can  transform  it  to  its  row  echelon  form 

fc  = Mb  = [°] 

by  a nonsir.gular  matrix  M.  (Here,  the  zero  is  an  (n  - l)-vector.)  Setting 

x = Mx 

we  get  a transformed  problem  analogous  to  (1 .!)-(! .2)  with 


J 
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■ — 


L L 


ft  = MAM'1  , "0  = (M1  )'1QM'1  > 0, 

= b'QB  > 0,  and  u = u. 

Let  us  assume  now  that  such  transformations  have  already  occurred  so  that 


(3.2) 


b = [°] 


and,  using  corresponding  partitioning,  write 


(3.3) 


A = 


An  ai2 


a21  a22 


, Q = 


Qn  qi2 


X = 


, and  p = 


12  q22 


where  a12,  a^,  q.^,  x-j  and  p^  are  (n  - l)-vectors  and  a22,  q22,  x2, 

and  p2  are  scalars  with  q22  > 0.  The  canonical  equations  then  take  the  form 


(3.4) 


X1  Allxl  + a12x2 


x2  = a^x-j  + a22x2  + u 


P1  ~ “Qllxl  " q12X2  Allpl  " a21p2 


and 


p2  _ql 2X1  “ q22x2  _ a12pl  " a22P2 


subject  to  the  end  conditions  that  x(0)  is  prescribed  and  p(T)  = 0. 


--- 
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With  the  control  constraint  |u|  <_  m,  the  maximum  principle  implies  that 
the  control  is  either  saturated  or  u = . Since  the  unconstrained  problem 

has  an  optimal  control  which  is  initially  unbounded  and  then  follows  a singular 
arc  as  6 ->  0,  it  is  natural  to  seek  a bounded  control  which  is  saturated  on  an 
initial  interval  0 <_  t <_  t-j  and  unsaturated  for  t-j  < t <_T.  For  such  switching 


solutions,  we 

must  satisfy  the  initial 

value  problem 

xi  = Anxi  + ai2x 

2 ’ x-j  (0) 

given 

(3.5) 

X2  = a21Xl  + a22X 

2 + m,  x2(0) 

given 

for  0 <_ 

t < 

t-j  and 

the  singularly  perturbed  two-poi 

nt  problem 

X1  = 

^llxl  + ai2X2’  x](t-j)  - 

x^tj) 

(3.6) 

= 6^a ' 
a21 

!X1  + £ a22x2  “ p2 

5 x2  ^1  ^ x2  ( ) 

pl  = 

-Qnxi 

- O-j 2X2  ‘ A-J  i Pq 

a21P2’  P1^T^ 

= 0 

P2  = 

-q|2»2 

- q22x2  “ al2pl 

a?2P2 * P2^1  ^ 

= +m£2,  p2(T)  = 0 

for  t-j  £ t £ T where  the  switching  time  tj  is  still  unspecified.  Besides 
requiring  continuity  of  the  states,  costates,  and  control  at  t-j , we  must  check 
the  control  remains  saturated  until  t-j  and  unsaturated  thereafter.  One  would 
expect  that  the  sign  of  the  initial  impulse  for  the  unconstrained  problem  would 
generally  predict  which  control  bound  would  initially  saturate. 

Our  previous  experience  with  singular  perturbation  problems  shows  that  the 
asymptotic  solution  on  tn  < t £ T will  be  of  the  form 
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x^t.S)  = X-j(t,£)  + t2rn-|  (t  ,t)  + fc2n-|(o,£) 

x?(t,£)  = X?(t,6)  + £m?(-r,£)  + £n„(0,£) 

(3.7)  2222 

P-jCt.e)  = P^t.6)  + ^(t.6)  + €2gi(a,€) 

p2(t,e)  = £2P2(t,e)  + e2f2(T3e)  + e2g 2(a,e) 

where  the  functions  of 

t = (t  - t^/e  or  o = (t  - t)/e 

tend  to  zero  as  that  "stretched  variable"  tends  to  infinity  (cf.  O'fialley  and 

p 

Jameson  (1975),  noting  that  we  need  p2  = 0(£  )).  Within  (t-|,T),  the  solution 
will  be  asymptotically  represented  by  the  outer  solution 

(3.8)  (XrX2,Pr£2P2) 

which  has  a power  series  expansion  in  t.  Its  leading  term  (xio’X20,P10’^  P20^ 
will  lie  along  a singular  arc  of  order  one.  We  note  that  the  initial  boundary 

2 

layer  term  f2(-r,0)  allows  the  nonuniform  convergence  of  the  control  u = -p2/£ 
as  £ 0 at  the  switching  point  t-|  from  +m  to  -P2Q(t^),  i.e.,  the 

jump  from  its  constrained  value  to  its  limiting  value  along  the  singular  arc. 
Analogously,  the  terminal  boundary  layer  term  g2(o,0)  allows  a jump  in 
p2/c  from  the  singular  arc  value  P2q(T  ^ t0  termina^  value  zero.  Explicit 
calculation  of  the  boundary  layer  terms  will  not  be  required. 

Since  (3.8)  must  satisfy  (3.6),  the  limiting  outer  solution  will  necessarily 
satisfy  the  reduced  problem 


15 


(3.9) 


fX10  A11X10  + a12X20’  XW  MV 
X20  = a21 X1 0 + a22X20  ' P20’  X20^V  = W 


P10  “Q1 1 X1 0 " q12X20  " A11P10’  PMT>  0 


0 = “q12X10  " q22X20  " al 2P1 0 


for  t-j  < t < T.  We  can  solve  immediately  for  X^g  and  P^g  as  functions  of 
X^g  and  P^g  (since  q^  i 0).  This  leaves  us  with  the  linear  two-point  problem 


(3.10) 


X10  ” AX10  " •SP1 0 ’ XW  = W 


^in  ~ '^in  ~ ^'^in’  ^in(^)  = 0 


i or  A - a]2q22q12’  ^11  - q12q22q12  — a^^  ^ a12q22a12  — 

As  usual,  it  is  convenient  to  solve  this  standard  regulator  problem  by  setting 


(3-11) 


P10(t)  = K(t)X10(t) 


where  the  (n  - 1)  x (n  - 1)  matrix  K satisfies  the  Riccati  differential 
equation 


(3.12) 


K + KA  + A'K  + d s KSK,  K( T)  = 0. 


Standard  arguments  (cf.,  e.g.,  Athans  and  Falb)  guarantee  the  existence  of  a 
unique  symmetric  solution  K >_  0 for  all  t < T.  Thus,  X^g  satisfies  the 
initial  value  problem 


(3.13) 


X1  g ~ (A  - SK)X^g,  X^g(t^)  X^(t-j) 


I 


while 
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(3-14)  X2Q  - ”^22^12  + a12^X10’  X20^1^  = x2^1^‘ 

The  initial  values  x^(t-j)  and  x ^ ( t -j ) must  be  determined  by  integrating  the 
saturated  control  problem  (3.5)  on  [0,t-,]  and  determining  the  switching  time 
t-j  by  solving  the  nonlinear  scalar  equation 

(3.15)  x2(t^)  = -q22(qj2  + a j2K(t1))x(t^). 


For  certain  simple  problems  (like  our  preceding  examples),  t-j  can  be  explicitly 
determined.  In  general,  however,  one  must  attempt  to  solve  (3.15)  for  its  least 
positive  root  t-|  numerically. 

To  check  that  our  solution  candidate  remains  saturated  on  0 <_  t < t-j , we 
integrate  the  linear  system 


f P1  = “Ailpl  ' a12p2  " Q11xl  _ qlZx2’  = P10^1'  = K(t|)xj(tj) 

(3.16)  ) 

Ip2  _ -ai 2P1  ” a22P2  ~ P12X1  ~ ^22x2 ’ ^2 ( ^"1  ^ = ^ 

backwards  from  tj , insisting  that  p2  i 0 on  0 < t < t^  (the  nonhomogeneous 
terms  are  known,  since  x^  and  x2  follow  from  (3.5)).  Likewise,  we  must  be 
sure  that  p2  remains  unsaturated  for  t > t-| , i.e.. 


(3.17) 


201 


I a21 X1 0 + a22X20  " X2(T 


< m 


for  t-j  < t < T. 

Our  procedure  then  shows  us  how  to  produce  candidates  for  the  optimal 
solution  of  Case  1 problems.  Examples  can  surely  be  found  (cf.  next  section) 
where  the  troublesome  possibilities  mentioned  above  eliminate  some  potential 
solution  candidates.  For  other  problems,  it  produces  suboptimal  controls  which 
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may  still  be  of  value.  Occasionally,  we  can  be  sure  that  a computed  solution 
is  optimal.  An  example,  on  an  infinite  time  interval,  is 

x.j  = + u,  x-|  (0)  = 0 

= -u,  x2(0)  = -0.5,  | u | 1 

J = if  x?dt  with  x-i(T)  = x9(T)  = 0 for  T free 
Ll  01  1 i 

\ 

(cf.  Johnson  and  Gibson  (1963)).  We  should  note  how  effective  our  method  is  in 
selecting  the  singular  arc  solution  and  that  it  relates  to  the  popular  technique 
of  synthesizing  a control  by  integrating  backwards  from  terminal  time. 

Algori thm 

We  have  used  the  above  asymptotic  analysis  to  construct  the  following  algo- 
rithm to  find  numerical  solutions  of  Case  1 problems. 

1.  Solve  the  Riccati  differential  equation  (3.12)  backwards  from  t = T 
to  t = 0.  [We  used  Gear's  code  (1971)  for  all  of  our  numerical  inte- 
grations; however,  any  other  good  code  would  have  sufficed.] 

2.  Decide  which  way  the  control  bound  will  initially  saturate  and  select 
an  initial  guess  for  the  switching  time  t-j . 

3.  Integrate  equations  (3.5)  from  t = 0 to  the  latest  guess  for  the  switching 
time.  We  repeat  this  step  until  equation  (3.15)  is  satisfied  to  a suf- 
ficient degree  of  accuracy.  We  generate  successive  guesses  for  t-| 

using  a Newton-like  procedure  due  to  Brown  (1969).  If  a negative  root 
is  found,  it  can  be  eliminated  by  deflation  so  that  a second  iteration 
can  be  attempted.  Also,  if  the  iteration  is  diverging  the  sign  of  the 
control  bound  can  be  reversed  and  the  procedure  repeated. 


H ..  ■ 


Once  the  switching  time  has  been  found,  we  integrate  equations 
(3.16)  backwards  from  t-|  to  determine  the  costate  vector  in  the 
saturated  region  [0,t-j].  We  then  integrate  equation  (3.13)  forward 
from  t-j  and  use  (3.14),  (3.11),  and  the  second  of  (3.9)  to  determine 
the  solution  in  the  unsaturated  region.  During  the  integration,  we 
check  that  the  appropriate  bounds  on  p2  are  not  exceeded. 

Finally,  calculate  the  cost  by  integrating  (1.2)  using  the  trapezoidal 
rule. 


4.  TWO  MORE  DIFFICULT  PROBLEMS 
a.  Sirisena  (1970)  considered  the  Case  1 problem 


(4.1) 


x i Xi (0)  2 

x2  = u,  x2(0)  = 0 

J = \ C(xi  + Xo ) ^ + £2u2]dt,  | u | < 1 

n 1 c 


for  € = 0.  If  we,  instead,  determine  the  asymptotic  solution  as  6-^0  as  in 
Section  3,  we  determine  the  switching  time 


t1  = S5  - 1 « 1.236  (for  T » t] ) 


where  the  control  changes  from  u = -1  to  the  singular  arc.  Specifically,  the 


limiting  control  is 


u = t-j  e 


(t-t, ) 


for  t-|  < t < T 


and  the  corresponding  limiting  cost  is 


J = ^5(25/5  - 41)  * 0.993447. 
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This  control  would  seem  to  compare  favorably  to  Sirisena's  bang-bang  solution 

for  the  fixed  endpoint  problem  with  x(T)  = 0,  achieved  with  an  infinite 

number  of  switchings  with  first  switching  point  1.227  and  cost  0.993455.  If 

we  take  our  constraint  seriously,  however,  we  must  rule  out  this  control  since 

u(t^)  = /5  - 1 exceeds  the  bound  u = 1,  (see  Figure  3). 

Our  nearly  good  solution  suggests  that  we  might  instead  seek  a solution 

k where  the  control  switches  from  u = -1  on  some  interval  0 <_  t < t^  to 

u = 1 on  t-j  < t < tg  after  which  u follows  a singular  arc.  An  unsaturated 

control  must  be  constructed  to  smoothly  switch  from  -1  to  1 near  t^  (i.e., 

2 2 

from  p^  = £ to  -£  ) . This  type  of  interior  (or  transition)  layer  differs 
somewhat  from  our  earlier  boundary  layers  because  the  desired  limits  occur  for 
finite  negative  and  positive  values,  -n  and  n+>  of  the  appropriate  stretched 

variable  n = (t  - t^ )/€  , not  for  n = Moreover,  the  limiting  solution  is 
unrelated  to  a singular  arc.  For  these  reasons,  the  local  behavior  is  determined 
by  a regular  perturbation  procedure  on  (-n  ,n.)>  although  the  solution  will 

~ T 

still  exhibit  singular  perturbation  features  (i.e.,  nonuniform  convergence) 

■ 

because  the  t interval  corresponding  to  (-n  ,n+)  vanishes  as  6 -*  0.  As 
before,  the  boundary  layer  at  t ^ 1S  studied  by  a singular  perturbation  analysis 
which  allows  the  control  to  switch  from  its  constraint  set  to  a singular  arc. 

For  (4.1),  the  canonical  equations  are 
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(xrx2,u)  = (2  - 

2 2 2 2 2 

For  n_  <.  t - t-j  <_  £ n+>  we  have  u = -p2/«  and  p2 ( +€ crv_)  = + £ , 

P2(t])  = 0.  Moreover,  we'd  locally  obtain  x1 , x?,  pp  and  p2  as  power 

series  in  £ with  x-j , x2>  and  p^  having  nearly  the  constant  values 

xl(t-|),  x2 ( t t ) » anc*  P2(t-j).  Since  the  limiting  value  of  p2  jumps  sym- 
metrically from  -n_  to  n+»  we  have  n_  = n+  = ( | x-j  (t^)  + x2(tj)  + p^tpj)"1, 
provided  the  denominator  is  nonzero.  Proceding  to  t|  < t < t2,  we  must  have 
the  limiting  solution 

(xrx2,u)  = (^1  ( t] ) + *2(tJ)(t  - tj)  + (t  - tf)2/2,  x2(tj)  + t - t|,l) 


_ O O 

and  p2(t2)  = £ = P2(t-|).  For  t > t2,  we'll  have  u = -p2/£  and  (4.2)  implies 

p 

that  the  limiting  solution  (xiq’X20’P10’^  p2q)  must  satisfy  the  limiting 
singular  arc  problem 


X10 

o 

C\J 

X 

II 

X • 
no 
o 

II 

o 

pio 

0 

X 

1 

II 

o 

II 

1 

X 

o 

1 

p 

• X 
1 

II 

o 

<NJ 

P10’ 

p10  ” 

X1 0^2  ^ = X1  ( t2) 


X20^2^  ~ x2^2^ 


Thus,  X20  = _X10  " P10’  P20  = ~X20  = _X10’  ancJ  t*iere  remains  the  two  point 

problem  X-jq  - -X-|q  ~ piq>  P10  ~ P10*  X10^2^  ~ x1^t2)»  p10^^  = lienee, 

we  follow  the  singular  arc 


(4.3) 


/ -(t-t2) 

X10(t)  = x1(t2)e 

( 

^P1q( fc)  = °’  X20  = P20  = 'X10 
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for  t > t^.  In  order  for  the  control  to  remain  unsaturated  there,  we  must  have 


I x1 ( t~)  | < 1. 


Lastly,  we  obtain  and  p2  for  t < t2  by  integrating  the  costate  equations 
of  (4.2)  backwards  with  the  limiting  boundary  values  P-j  ( t2 ) = 0,  P2(t2)  = 0. 
Since  x-j  and  x2  are  defined  differently  for  0 <_  t < t-|  and  t-j  < t < t2, 
we  also  have  different  expressions  for  p-j  and  p2  on  these  intervals.  Our 
solution  is  now  completely  determined  up  to  specification  of  the  switching  times 
t-|  and  t2,  0 < t-|  < t2- 

The  switching  times  follow  from  the  two  scalar  equations 


(4.4) 


which  reduce  to 


*2^2)  ~ ^2^2^ 
and 

P2  ^ ^"i ) = 0 

l 


^<At)2  + (1  - t-j ) At  + (2  - t1  - [t2)  = 0 

3 pH 

l(At)3  + 1(3  - 2t-j ) (At)2  + ^{3  - 2t1  - lt2)At  + (2  - t1  - It2)  = 0 

for  positive  At  = t2  - t^.  Solving  these  equations  simultaneously  we  obtain 
t-j  = 1.227  and  t2  = 1.564  and  a resulting  cost  0.993455.  We  further  check 
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that  has  no  root  other  than  t^  on  0 < t < t2  and  that  | p2 ( t^ ) | 

= | x-j ( t^)  I <1-  We  note  that  these  results  (shown  in  Figure  4)  are  in 
agreement  with  those  of  Sirisena  (1970),  who  needed  an  infinite  number  of 
switchings.  Presumably,  an  analogous  procedure  could  be  used  for  other  Case  1 
problems  requiring  a finite  number  of  switchings. 

b.  Jacobson  et  al«  (1970)  considered  the  Case  2 problem 


(4.5) 


fxl  x2  ’ 


x2  = u. 


If5, ..2 


2.2, 


J(e)  =2 j (xi  + e u )dt. 


x-|  (0)  = 0 
x2(0)  = 1 
!u|  < i 


with  € = 0.  Since  the  simplest  Case  1 solutions  involve  a control  which 
switches  from  being  saturated  to  lying  on  a singular  arc,  we  might  attempt  to  find 
Case  2 solutions  which  switch  once  between  the  two  control  bounds  and  then  follow 
the  singular  arc.  (For  this  problem,  it  will  be  impossible  to  find  a control 
which  switches  directly  from  its  bound  to  the  singular  arc.)  Noting  that  the 
unconstrained  problem  has  a large  negative  initial  control  as  £ -*  0 (cf. 

O'Malley  and  Jameson  (1976)),  we  seek  a solution  such  that  u = -1  for  0 <_  t 

p 

< t-j , u = 1 for  t-j  < t < t2>  and  u = -p2/€  for  t near  t-j  and  for 
t2  < t < 5.  On  (t2,5),  the  limiting  solution  (x10’X2o’pio,p20^  wil1 
satisfy  the  limiting  canonical  system 


(4.6) 


x10  x20*  0 _p20’  p10  = ~x10’  and  p20  " ~p10’ 


! 


i.c.>  it  will  follow  the  trivial  second  order  singular  arc.  Integrating  the 
state  equations  (4.5),  we  obtain 
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(x1>x2^  = - ^t^,l  - t)  for  0 < t < t-j 

and 

(X1  >x2)  - ( x i ( 1 1 ) + x2 ( 1 1 ) (t-ti ) + 2*C "t — 1 1 ) , x2(t-| ) + ( t — 1 1 ) ) for  t-j  < t < t2- 

As  for  the  preceding  problem  the  transition  layer  at  t-|  and  the  boundary  layer 
at  t2  will  allow  discontinuities  in  u,  but  not  in  the  states  x-j  and  x^. 

Thus  we  must  have  x^  (t^)  = x^  (t^)  and  x^(t~)  = x^U^)  = 0,  i = 1,  2. 

These  conditions  imply  that 

(4.7)  t]  = 1 + — - and  t = 1 + J2. 

The  resulting  solution  (see  Figure  5)  compares  very  favorably  to  that  of  Jacobson 
et  al  ; the  cost  being  0.269  compared  to  their  value  0.277. 

Regrettably,  however,  the  control  just  computed  is  suboptimal . It  does  remain 
unsaturated  on  (t2,5) , but  the  costate  vector  must  satisfy  the  limiting  terminal 
value  problem 


( 


-Xp  p^t')  = 0 


(4.8)  / 

p2  = "Pi  > p2(*2)  = 0 

i 

for  t < t2-  Since  x-|  > 0 within  (0,t2),  p-j  and  p2  are  also  positive 
there,  so  we  can't  achieve  the  optimality  condition  p2(t-|)  = 0.  (Our  algorithm 
yields  p2(t|)  = 0.0104).  This  closure  problem  occurs  because  our  switching 
requirements  completely  determine  the  state,  costate,  and  switching  times.  The 
additional  optimality  condition  is  an  extra  constraint  without  a corresponding 
degree  of  freedom. 
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It  should  be  realized  that  the  optimal  control  for  this  problem  presumably 
involves  an  infinite  number  of  switchings  which  accumulate  at  a point  after 
which  the  control  follows  a singular  arc  (cf.  Fuller  (1963),  Grensted  and 
Fuller  (1965),  and  Marchal  (1973)),  although  fewer  switchings  may  be  possible 
for  special  initial  values.  Our  suboptimal  results  seem  very  good,  however. 

They  relate  to  Sirisena's  concept  of  an  optimal  r-switch  solution  (cf.  Sirisena 
(1970)). 

A related  Case  2 example  with  a trivial  singular  arc  requiring  an  infinite 
number  of  switchings  is 


x1  = x2. 


( x?  = u, 

r4 


J = 


1 


2 


2,^ 

x-|dt. 


x-,  (0)  = 2 
x2(0)  = -2 
|u|  < 1 


(cf.  Marchal  (1973)).  When  we  try  to  calculate  a one  switch  solution  starting 
at  u = -1 , we  get  t-|  = 2,  At  = t2  - t-|  = 0,  and  cost  1.6.  The  optimal 
solution  has  u = 1 until  t = 0.057,  switches  to  u = -1  until  about  t = 2.6, 
and  then  has  an  infinite  number  of  switchings  accumulating  at  t = 3,43  after 
which  u = 0.  The  optimal  cost  is  1.52.  Our  suboptimal  solution  is  presumably 
not  as  good  as  for  the  preceding  problem  where  we  were  able  to  switch  once  before 
meeting  the  singular  arc. 


5.  COMMENTS  ON  HIGHER  ORDER  SINGULAR  ARCS 
Experts  (cf.  e.g.,  Krener  (1976))  feel  that  infinite  switching  is  generic 
for  problems  with  singular  arcs  of  order  p,  P >_  2,  although  finite  switching 
may  be  possible  for  special  initial  conditions.  Our  preceding  discussion  supports 
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this  opinion.  For  such  problems,  we  can  readily  identify  the  appropriate  form  of 
singular  arc  solution  through  our  singular  perturbation  analysis,  and  can 
generally  find  a suboptimal  p - 1 switch  solution  converging  to  the  singular 
arc.  Suboptimal  solutions  with  more  switchings  (a  la  Sirisena)  might  also  be 
sought. 

Specifically,  for  Case  2,  we  have 

b'Qb  = 0 and  b'A'QAb  > 0. 

The  problem  can  be  transformed  so  that 


( •>  r 


o' 

A11  a 12  ° 

_a 

— i 

ro 

o 

(5.1)  b = 

0 

, A = 

a21  a22  1 

, and  Q = 

q12  q22  0 

1 

a31  a32  a33 

0 0 0 

where  A-j-j  and  are  matrices  of  dimension  (n  - 2)  x (n  - 2)  and  q^2  1S  a 

positive  scalar  (cf.  O'Malley  and  Jameson  (1976)).  Introducing  corresponding 
state  and  costate  vectors 


X1 

P1 

X = 

x2 

and  p = 

P2 

X, 

l 3J 

P3. 

o 

the  scalar  optimal  control  either  lies  on  its  bound  |u|  = m or  else  u = -p3/£  . 
On  the  singular  arc,  the  control  should  be  unsaturated  and  the  canonical 
system  has  the  form 


and  there  remains  the  linear  2(n  - 2)  -dimensional  problem 


(X10  (All  “ a12q22q12^X10  ' a12q22a12P10 
(5.5)  j 

P10  = "(Q11  “ q12q22q12^X10  " ^All  ' q12q22a12^P10 

analogous  to  (3.10).  Since  Q >_  0 and  q^  > 0 imply  that  - 1]2q22q12  — 0> 

standard  results  imply  that  the  two  point  problem  with  X,q(1:)  prescribed  and 
p1q(T)  = 0 will  have  a unique  solution  on  any  interval  [^,T].  Thus,  the 
singular  arc  solution  with  control  u -P^q  can  be  easily  computed  on  any 
such  interval  provided  the  state  x,(t  ) is  known  from  the  initial  trajectory 
[o.t)  with  bang-bang  control.  In  particular.  Case  2 problems  converge  to 
the  trivial  singular  arc  solution  when  n = 2.  The  Case  2 problem 
• • • 

x-j  — *2  ~ x3’  x3  ~~  ^ ^ ^ ( 0 ) — I*  x^lO)  ■“  0,  ( 0 ) “ 0 

|j  = l|0(xi  + x2^dt>  iul  i 1 

will,  for  pvample,  have  a nontrivial  singular  arc  (cf.  Sirisena  (1970)). 


We  car-  analogously  determine  the  singular  arcs  appropriate  for  Case  p 
problems,  p > 2.  For  example,  we  find  the  trivial  singular  arc  for  the  Case 
3 problem 
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x -j  ~ X2 j x 2 x^>  x3 


x^dt,  | u | <_  1 


u,  x-j  (0)  = 1,  x2(0)  = 0,  x3(0)  = 0 


(cf.  Grensted  and  Fuller  (1965)  and  Sirisena  (1974)).  The  optimal  solution 
switches  infinitely  often  before  landing  on  the  singular  arc  and  has  an  optimum 
cost  1.2521.  Our  simplest  suboptimal  solution  switches  at  t-(  ~ .794  and 
3t-|  and  reaches  the  singular  arc  at  4t^  with  a cost  of  1 .2665. 
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